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INTRODUCTION 


The  classical  bound  on  the  error  in  linear  interpolation  of  function  f  on 
interval  (a,b)  is  given  by 

where 

ilf"ll(a,b)  =  sup  i  f’’(x)i 
a<x^b 

We  wish  to  obtain  meshes  (x^ ,  1  <  i  $  n)  whicn  will  somehow  equilibrate  the 
"error"  over  each  subinterval. 


GOOD  MESHES 

C.  deBoor  (ref  1)  has  supplied  us  with  a  computationally  simple  method  for 
generating  what  he  calls  "good"  meshes.  His  idea  is  to  make  the  classical  bound 
roughly  constant 

|(Xi+i-Xi)*llf"II(x.j,Xi-+i)  =  constant  1  ^  i  <  n 
This  is  equivalent  to 

(Xi+i-XT)llf"ll(x.i,Xi+i)  =  constant 


or 


.xi+i) 


dx 


c 


As  n  becomes  large  and  x.j+i-xi  -*  0  for  all  i,  we  approximate  the  integrand  by 
I  f"(x)|  ^  =  g(x). 

Thus,  we  have  a  relatively  simple  problem  to  solve 


’‘i  +  l 

f  g{x)dx  =  c  1  $  i  <  n 
Xi 

/X 

If  we  define  G(x)  =  /  g(t)dt,  we  have  G(x^-)  =  (i-l)c,  therefore 


1 


G(Xi) 

G(x“] 


--7  and  xj 
n-1  1 


G“’ 


1-1 

n-i 


S{xn) 


1  < 


<  n 


In  practice,  we  may  only  have  a  positive,  continuous,  piecewise  linear  estimate 
of  g  over  some  mesh  u.  We  denote  this  estimate  of  g  by  v.  G  as  defined  by 

G(x)  =  /  v(t)dt  <  X  <  Urn 

would  then  be  piecewise  quadratic  and  invertible  in  the  following  manner: 


G-’ (G*) 


X* 


u,  +  2 


(G*-Gt) 

v^+yo 


where 


G^  =  0,  Gj+i  =  Gj  +  (uj+i-uj) (Vj+Vj+i)/2  1  ^  j  <  m 


G  -j  ^  G*  ^  G  j 
p  =  (G*-Gi)/(Gi+i-Gi) 
and 

2  2 

D  =  (l-p)v^  +  pvi+i 

All  this  is  common  knowledge.  Unfortunately,  good  meshes  do  not  always  seem 
quite  as  good  as  we  might  like  them  to  be.  The  shorter  subintervals  have  a 
fairly  uniform  error  bound  pattern,  but  the  lengths  of  the  longer  subintervals 
are  overestimated,  yielding  larger  error  bounds.  This  is  due  to  the  fact  that 
the  integral  of  the  norm  of  g  is  underestimated  by  the  integral  of  g.  In  fact, 
it  is  easy  to  prove  that  for  f(x)  =  x^(p  >  2,  0  ^  x  ^  1),  the  largest  error 
bound  on  a  good  mesh  is  exactly  equal  to  the  largest  error  bound  on  the 
corresponding  uniform  mesh  (xi+j-x-j  =  const). 

In  order  to  get  what  we  might  call  "better”  meshes,  we  go  back  to  the 
original  problem: 
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Find  n-2  x's  (xj  and  x^,  fixed)  such  that 


(Xi+i-Xi)llgll(x^,XT  +  i)  =  "9n{xi,Xi  +  i)dx  =  c  1  $  i  <  n 

This  problem  is  described  in  Reference  1  as  being  rather  difficult  to  solve  in 
general.  Even  if  we  knew  what  c  was,  solving 

(xi+i-x-i )  llgll  ,x.j+i )  “  c 

for  x^+i  given  x^  would  not  be  easy.  The  problem  obtained  upon  substituting  v 
for  g,  however,  is  quite  easy  to  solve  (ref  2).  In  addition,  if  v  is  a  very 
good  approximation  to  g  (with  m  >>  n),  we  get  a  virtually  constant  error  bound 
for  the  entire  mesh! 

We  refer  to  the  following  equation  as  the  "stepping"  equation: 

(/3-a)  llvll  (a,j3)  =  c 

The  solution  of  the  stepping  equation  for  /3  given  a  and  c  represents  the  central 
part  of  our  algorithm  for  obtaining  better  meshes.  Although  the  stepping 
equation  is  nonlinear  in  /S,  the  piecewise  linearity  of  v  enables  us  to  solve  it 
noniteratively .  For  given  a  and  c,  we  solve  the  stepping  equation  in  the 
following  manner.  Suppose  a  e  (u^.u-j+j)  and  we  have  located  j  such  that 

(Uj-a)llvll(a,uj)  <  c  <  (Uj+l'®)  "vH  (a,  uj+j ) 

Hence, 


0  €  (Uj,Uj+i) 

To  locate  j,  we  simply  search  from  left  to  right,  computing  the  norms  as  we  go 
and  checking  the  previous  inequality.  We  use 

"vll(a,u.j+l)  =  Max{v(a)  ,Vt+i} 

and 


llvll 


(a,Uj+i) 


Max(llvll(a^u^)  ,vj+i) 
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Now  if 

llvll(a,uj)  =  llvll(a,uj  +  i) 

we  have 

(;3-a)llvll(a,Uj)  =  c 

hence 

^  =  a  +  c/llvll{a^u^) 

But  if 

"'/"(a.Uj)  <  (o,Uj+i) 

there  is  a 

t  e  (Uj,uj+i) 

such  that 

'•'^"(a.uj)  =  llvllja^t)  =  v(t) 

and  this  t  is  given  by 

t  =  uj  +  (llvll(a,uj)-vj)/s 

where 

S  =  (Vj+i-Vj)/(Uj+l-Uj) 

Now  i  f 

(t-a)llvll(a^uj)  >  c 
must  lie  to  the  left  of  t  and 

0  =  a  +  c/llvll(a,u^j 

as  before,  but  if 

{t-a)llvll(a^Uj)  <  c 


/3  lies  to  the  right  of  t  and 

{P-a)llvil(<j^0)  =  c 

However,  in  this  case, 

ilvll(o,^)  =  v(/3)  =  vj  +  s()3-Uj) 
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Therefore , 


(0-0) (Vj+S(0-Uj) )  =  c 
or 

(0-0)  (vj+s(0-a- o-'jj) )  =  c 
or 

s(0-o)*  +  k(0-o)  -  c  =  0 


where 


k  =  Vj  +  s(o-uj) 

Solving  this  simple  quadratic  equation  for  0  -  o  yields 

0=0+  (/k2+4sc-k)/(2s)  for  k  <  0 

and 

0=0+  2c/(k+/k*+4sc)  for  k  >  0 

Having  elaborated  the  solution  to  the  stepping  equation  for  arbitrary  c,  we 
now  consider  obtaining  the  correct  value  of  c  by  defining  the  function  u 

p(c)  =  y-n 

where  v  is  the  number  of  x's  we  get  by  solving  the  stepping  equation  y-2  times 
(xj  and  Xy  =  Xp  being  fixed).  It  is  intuitively  clear  that  for  small  c,  p  will 
be  positive  and  for  large  c,  p  will  be  negative.  Since  u  is  a  step  function,  we 
are  interested  only  in  its  leftmost  zero,  the  correct  valuf  of  c.  When  the 
correct  value  of  c  has  been  obtained,  we  shall  have  concurrently  obtained  the 
better  or  uniform  error  bound  mesh.  Therefore,  we  see  that  to  find  the  better 
mesh,  we  need  only  solve  a  single  nonlinear  equation  in  a  single  unknown  (by 
modified  bisection),  where  each  evaluation  of  p  involves  an  0(m)  search  through 
the  (u,v)  data  and  0(n)  solutions  of  simple  linear  or  quadratic  equations. 
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Suppose  we  are  not  satisfied  with  better  meshes  and  decide  to  take  the 
additional  step  of  finoing  the  best  possible  mesh.  We  define  this  brst  mesh  as 
one  in  which  the  maximum  absolute  error  is  constant,  independent  of  subinterval. 
Now  we  must  work  in  terms  of  exact  error  instead  of  error  bounds.  The  exact 
error  in  linear  interpolation  of  f  on  interval  (a,b)  is  given  by 

e{x)  =  (t-a)f"(t)dt 

D-a  ^ 


*  (b-t)f"(t)dt 

□  -a 


If  x^  is  the  maximizing  or  minimizing  point  of  e, 

e' (Xm)  =  0 


implies 

/  {t-a)f"(t)dt  =  ;  (b-t)f’'(t)dt 

a 

which,  in  turn,  implies  that 

e(^m)  =  (t-a)f"(t)dt  =  'b-t ) f " { t )at 
a  ^m 

These  last  two  equations  tell  us  that  for  given  a,  E,  and  f",  we  can,  in  prin¬ 
ciple,  solve 

E  =  /  *"  (t-a)f"(t)dt 
a 


for  x^  and  then  solve 

E  =  /^  (b-t)f’*(t)dt 

for  b.  We  therefore  see  that  finding  the  best  mesh  is  not  a  very  different 
process  from  finding  the  better  mesh.  The  major  difference  is  chat  two 
equations  must  be  solved  in  the  stepping  process  instead  of  one.  This  stepping 
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process  is  further  simplified  to  solving  simple  quedratic  equations  if  we  use  a 
piecewise  constant  aoproximation  to  f". 

The  first  half  of  the  stepping  process  amounts  to  solving 

R(j9)  =  +  E  for  /3 

g-.'en  a  where 

jS 

R(/3)  =  /  (t-a)f”(t)dt 
a 

if  f"(x)  =  on  (u-ijU^^^)'  have 
,  ,  ^i 

R((3)  =  /  (t-a)f’'lt)dt  /  (t-alc^dt  (iSe  ( u^ ,  u^ +  -,  ) ) 
a 

P,j3)  =  R(u^)  +  ^  CT(t-a)2  !  ^ 

=  R(Ut)  +  I  Ct  (  (/3-a)  ^-(u^-a)  2 ) 

=  R(u^)  +  I  Ci(/3-u-j){/3+UT-2a) 

The  recursion  ^or  R  is 

^i+1  =  f^i  ^  Ci  (Ui+^-Ui )  (uT+u,  +  i-2a) 

The  second  half  of  the  stepping  process  amovnts  to  solving 

S(y)  =  +  E  for  y 

given  where 

S(y)  =  (y-t)f"(t)dt 

0 

For  y  e  (Uj,u^^^),  we  have 

U-j  y 

S(y)  =  /  (y-t)f"(t)dt  +  /  (y-t)c-iat 
a  u^ 

*1  • 

=  /  (u^  t+y-ui)f"(t)dt  -  I  CT(y-t)^l  ^ 

=  S(u^)  +  (y-ui)  /  f”(t)dt  +  3  CT(y-ui)2 

P  ‘ 
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Letting 


J  f(t)dt 
/3 

we  have 

fi+l  =  Ti  +  CT(Ui+l-Ui) 

and  the  recursion  for  S 

^i+1  “  S-j  +  (u-j  +  2~^i)^i  2  ^  i  ( ^i +  1“*^  i )  ^ 

We  mention  in  passing  that  after  we  have  computed  j3,  and  as  we  search  for  y,  we 
will  need  to  compute  additional  /3's  if  there  are  inflection  points  present  in  f. 

Having  given  this  rather  brief  sketch  of  the  best  mesh  process,  we  proceed 
to  some  interesting  graphical  results  shown  in  Figures  1  through  4.  We  have 
obtained  good,  better,  and  best  meshes  for  the  simple  test  functions  x*(l-x)* 
and  X’*,  and  we  may  describe  our  little  odyssey  as  follows.  Starting  with 
deBoor's  good  mesh  with  its  predictably  large  bounds  on  the  long  subintervals, 
we  then  proceed  to  flatten  these  bounds  out  almost  perfectly  through  the  better 
mesh.  Proceeding  one  additional  step  to  obtain  the  best  mesh  (whose  maximum 
error  can  indeed  be  seen  to  be  constant),  we  then  note  the  striking  similarity 
oetween  the  good  mesh  and  the  best  mesh.  One  might  therefore  say  that  in  going 
from  good  to  better  to  best  meshes,  we  have  come  nearly  full  circle  and  can  see 
the  wisdom  of  the  old  French  proverb;  ’’The  more  things  change,  the  more  they 
remain  the  same."  It  would  seem  that  the  only  thing  deBoor  can  be  criticized 
for  is  excessive  modesty  in  referring  to  his  meshes  as  merely  good  when,  in 
fact,  they  are  nearly  best. 
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Figure  3.  Second  test  function  defined  on  good, 
better,  and  best  meshes. 
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